We use RStudio IDE. Create a .R file (ours is called webinar1_2020_07_29.R) and write the following lines of code. Ctrl-Return to run the line that cursor is on or run all highlighted lines. Or use the button in top right part of the editor.
To add functionalities to your .R script you must load the specific libraries. REMEMBER: make sure that they are installed in your system. If not you can install through RStudio on menu => Tools => Install Packages or directly through R console with command {r eval=F } install.packages("NAME OF PACKAGE")
Below we add three libraries which add functionalities that we need
My custom Coordinate Reference System (CRS) projection, Lambert Conical Conformal (LCC), NOT secant but tangent at lat=45.827 and lon=11.625 - for more info:
My points (in lat long) over the regular grid/lattice these points are 666.67 m apart
## Reading layer `points' from data source `/archivio/R/shared/webinar/data/webinar1_2020_07_29/points.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1315 features and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11.83267 ymin: 46.33037 xmax: 12.14322 ymax: 46.55611
## geographic CRS: WGS 84
Convert my points from a Geographic (latitude and longitude degrees) to my custom CRS projection
Points to a square Polygon area. How? It is trivial, but faster way was to create a rhomboid with a certain distance from the point (“radius”) using a buffer of half the distance between points (333.33 m) . The minimum bounding box of the buffer with that radius will be a square with sides 2x the size of the radius of the buffer.
NB: nQuadSegs is 1 which means a chord(segment) at each quadrand, i.e. 1/4th of a circunference; this give a rhomboid which is a “very very rough circle” - it saves a lot of memory from drawing with a defaul value of nQuadSegs=30 - it makes a difference when processing a large number of points.
We buffer each feature (point) and view results
Then we apply a function to each “rhomboid” for creating the square. It is a bit more complex as the function includes another function inside.
st_bbox_by_feature = function(geom) {
## make sure that object "geom" becomes a geometry object;
## geom is your data set with all the single geometries/polygons
geom2 = st_geometry(geom)
#' Function to:
#' (a) take a single geometry,
#' (b) create a bounding box with st_bbox and
#' (c) convert the bbox object to a new geometry (a square)
f <- function(single.geom) {
st_as_sfc( st_bbox(single.geom,), crs=myproj)
}
#' This line calls the function above LOOPING over each single geometry (lapply)
#' in the geom2 object
do.call("c", lapply(geom2, f))
}This line calls the above function over our data (rhomboids)
Assign the CRS to the new dataset as it got lost. To check the CRS of our data we can use {r} st_crs(tiles)
We view all results all over the same map (small subset of data - the first 20 features).
Convert to Latitude and Longitude and save result to shapefile for future upload to Google Earth Engine
Import shapefile files “tiles.XXX” by left panel “assets” and click “NEW” and choose shapefile from dropdown menu.
DO THINGS in GEE
code shared https://earthengine.googlesource.com/users/2020_Kanan/summer_webinar_series in script https://earthengine.googlesource.com/users/2020_Kanan/summer_webinar_series/+/636ab64efd7b09a0d45d3d10a992b8a44bd119c7/webinar1
EXPORT THE RESULTS TO A SHAPEFILE CALLED tiles_withData.shp
Read the file exported from GEE
## Reading layer `tiles_withData' from data source `/archivio/R/shared/webinar/data/webinar1_2020_07_29/tiles_withData.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1315 features and 6 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 11.82833 ymin: 46.32735 xmax: 12.14759 ymax: 46.55912
## geographic CRS: WGS 84
Create a new attribute column with the interquartile range (The difference between p75 and p25, i.e. 75th and 25th percentile)